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Abstract 

We develop a variational method to obtain accurate bounds for the eigenenergies of H = — A + V 
in arbitrary dimensions N > 1, where V(r) is the nonpolynomial oscillator potential V(r) = 
r 2 + ^ 2 ; ^ € ( — 00 i 00 )j3 > 0. The variational bounds are compared with results previously 
obtained in the literature. An infinite set of exact solutions is also obtained and used as a source of 
comparison eigenvalues. 
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1. Introduction 

The following is a variational study of the eigenvalues of the Hamiltonian H = A + V(r), where 
the non-polynomial oscillator potential V(r) is given in arbitrary dimensions by 

V(r)=r 2 + T ^, 5 >0,Ae(-oo,oo). (1.1) 

Various techniques may be found in the literature [1-53] for the estimation of the spectrum corre- 
sponding to this potential. Some of these techniques have been developed only for special values of 
A and g . The purpose of the present work is to introduce a unified variational technique to compute 
the eigenvalues of Schrodinger's equation with potential (1.1) in arbitrary dimension N > 1 and for 

arbitrary angular quantum number I = 0, 1, 2, The variational method presented here has the 

advantage of being valid for arbitrary values of A , g > in arbitrary dimensional N . The article is 
organized as follows. In section 2, we introduce a variational technique based on the Gol'dman and 
Krivchenkov Hamiltonian as a solvable model. Thereafter, we use its eigensolutions to compute the 
matrix elements for the operators r 2 /(l + gr 2 ) , g > . Closed- form analytic expressions in terms 
of the incomplete gamma function are obtained for the matrix elements. In section 3, application to 
the non-polynomial oscillator potential is discussed and we compare our results with those by other 
techniques avaliable in the literature. Some remarks regarding the perturbation expansions are also 
given in section 4. Our conclusions are presented in section 5. For comparison, we review, in the 
appendix, the exact solutions of the Schrodinger equation with the potential (1.1) valid for certain 
relations between the parameters A and g . These exact solutions are valid for A < and can be 
seen as generalizations of the exact solutions in the one-dimensional case developed earlier in the 
literature. The purpose of exploring these exact solutions is to use them for testing the accuracy of 
the variational method we develop in section 2. 

2. Variational approximation by means of an exactly solvable model 

The radial Schrodinger equation in N -space dimension for the n th quantum state of the non- 
polynomial oscillator potential (1.1) can be written, with %/2m = 1 , as 

h«w - ( - £ + *> + <LtMzim±Mzm + T ^)« M . « M , 

where I = for the N = 1 case. Here the parameter B is introduced to simplify the computation of 
the matrix elements of H ; it may be set equal to 1 afterwards, without disturbing the computation. 
For N > 1 , n refers to the radial quantum number (that is to say, n — 1 is the number of radial 
nodes). For certain relations between the parameters A and g, exact solutions can be obtained 
[22,23,31,32,35,36,38,39]. These exact solutions are summarized in the appendix and are valid for 
arbitrary values of 7 = I + ^ (N — 3), i.e. not restricted to certain integer values of N and I as 
reported in the literature. They will be used to verify the numerical computations by the variational 
technique we are about to develop. Our variational approach is based on the introduction of a 
dummy variational parameter A on which the full perturbation expansion does not depend (though 
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the truncated expansion does). The optimal value of the variational parameter is then determined 
by a minimization process. Let us write the Hamiltonian H in (2.1) as 

H = ffo + ( TT ^-^) I A>_-- (2.2) 

where 

g .--^ + ^ + g±MzMjJt^»±i ,,- M ,t,.. ,, 3) 

The unperturbed Hamiltonian H n , known as the Gol'dman and Krivchenkov Hamiltonian, is one 
of the few systems that admit exact analytical solutions [54] . It is the generalization of the familiar 
harmonic oscillator in 3-dimension —d 2 /dr 2 + Br 2 +1(1 + l)/r 2 , where the generalization lies in the 
parameter A ranging over [— j,oo) instead of only values determined by the angular momentum 
quantum numbers I = 0, 1,2, . . . . The energy spectrum of the Schrodingcr equation Hqi/j^ = e^i^nl 
is given, in terms of parameters A and B , by 



e^ = 2/3(2n + C), n = 0, 1, 2, . . . , (/? = Vb, < = 1 + yj A + (I + | - l) 2 
where the normalized wavefunctions are given explicitly by 



(2.4) 



f3r 2 , (n, 1 = 0,1,2,..., N > 1). (2.5) 



The confluent hypergeometric function 1F1 is defined by means of the polynomial 

1F1 ^ ™ z^j = - , ( n -degree polynomial in z ) (2.6) 

where the shifted factorial (a)„ is defined by 

(a) = 1, (a) n =a(a + l)(a + 2)...(a + n- 1), for n = 1,2,3, .. . (2.7) 

The singular basis, consisting of the set of exact solutions of Ho , serves as a good starting point for 
the analysis of the non-polynomial oscillator Hamiltonian H given in (2.1). Indeed, the use of the 
basis (2.5) will allow us to overcome the difficulties which arise if one uses Hermite functions [41]. 
The matrix elements of the Hamiltonian H , (2.1), in terms of the orthonormal basis (2.5), read 



2 



n\ry. = 2( In }6 mn + <^lj^l^> - <*&|A" 2 |,$>, (2-8) 



where the Kronecker delta d mn = for m ^ n and 1 for m = n . Consequently, to make use 
of the matrix elements (2.8), we need to evaluate the expectation values <V , mzl i+gr 2 IV^;> anc ^ 
< i'mi\ r ~ 2 \' l Pni > ^ means of the orthonormal basis (2.5). This is sunnarized by the following two 
theorems. 

Theorem 1: By means of the orthonormal basis {tpnl}?^=o > where ^(r) given by (2.5), we have 



<V m i\ 1 + gr2 ^m>-J l+gr2 dr-( 1) 5 c+ir(C)V n!m! 
o 

n 171 f ^ ( \ / q \ 

x 
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Here T{—Q — i — j;^-) is the incomplete gamma function defined by T(a,z) = J z t a 1 e t dt. 



Proof: Direct integration of J °° r ^"^^"'^ rfr using (2.5) yields 



«W rff = 2(-l)™+^ /(C)„(C), 



1 + gr 2 



r(C) 



n\m\ 



1 + <yr 2 



-iFi 



c 



(3r 2 1F1 



/3r 2 dr 



r(C) 



2(-ir+»/?< (C)n(C) m f f (-«)i(-m)# +j fr 2(+2i + 2 j+ l 



EE 



»!»'! (Oi(Oi<yi j i •/■- 



y i 



c 



(2.10) 

where we have used the terminated scries representation of the confluent hypergeometric function 
(2.6). The question is now reduced to the computation of the finite integral 



r 2C+2i+2j+l 

1 + gr 2 



r^dr = 4 



i r x i+i+j 



2 J 1 + gx 



~ 0x dx, (after substituting r 2 = x). (2.11) 



Such integrals have been questioned in the literature [5,9,14,28,41]; therefore, we shall give a detail 
calculation of it here. Differentiation of (2.11) with respect to (3 yields 



dl(l3) 1 f x C+i+j+ \-p Xdx 



(2.12) 

dp 2 J 1+gx y ' 

o 

Multiplying (2.12) by g and subtracting it from (2.11) implies 

oo 

m-g d -^ = \jx^e-^dx 



The integral on the right hand side can be easily calculated by using the definition of the gamma 
function T(z) — e~H z ~ 1 dt, Re(z) > , and yields 

This is an elementary first-order differential equation which can be easily solved as 



e » 



2 5 Je>0 



where the constant C is chosen such that I{j3) — > as (3 — > oo (a conclusion which follows from the 
definition of the Laplace transform). Therefore, we have, after simple algebra and use of L'Hospital's 
rule, that 



r(c + i + j + i) 



[ X e-h- {C+t+1+1) dv = r(C + 2 + J+ / ? ) r e-^+^+V 
h ' 2qt+ l +3+ l e-^ H 



dr\ 



2ge~9 J0 2gS 
which is the definition of the incomplete gamma function. This completes the proof. □ 
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Theorem 2: By means of (2.5) , we have 



<< l \r- 2 \^ l >={ 



C-i 



if n = ra, 



(2.13) 



Proof: Direct integration of J Q r 2 'ipni( r )' l l } ni( r )d r using (2.5) yields 



<^Z(r)\r- 2 \^Z(r)> 



2(-l)™+"^ /(C)„(C) 



^ ^ r^e-^^i-n; C; Pr^F^-m; C; (ir 2 )dr. 
l - Jo 



T(C) V nm\ 

The finite integral on the right-hand side follows directly by means of the identity [54] 



oo 



X-a-l e -Pr* iFi 



C 



(3r 2 1F1 



c 



2/? C -^(C)n 



(2.14) 



The proof of the theorem then follows from the fact that the hypergcometric function 

{-m) k {a)k{b) k 



. -m,a,b 
c,d 



1 =£ 



k=0 



(c)k(d)k k\ 



(m — degree polynomial). 



collapses to unity when a = 2 . □ 

Theorems 1 and 2 provide closed-form analytic expressions for the matrix elements (2.8); these allow 
us to obtain accurate upper bounds for the Hamiltonian (2.1) by means of the following optimization 
problem: 



Eni < min diag 



Hoo Hoi 



Hod-i 

HiD-1 



H mn = <ilCimZ>, (2-15) 



^-ff_D-10 Hd-u ■■■ Hd-\D-1- 

for fixed B . 

3. Applications and Numerical Results 

In order to avoid the numerical complication [55] of computing the incomplete gamma functions in 
(2.8), we have, by means of the recurrence relation [56] 



T(a + l,z) = aT(a, z) + z a e" 



that 



Thus 



r,, \ T(a + n,z) 
™ = —^ n 6 



E 



^a-\-k— 1 

\ (a)fc 



v( , ■ ■ f>, (-i) i+J r(-C,f) 



g\(+i+3 + l -S. 



1+3 



(-) 



e s 



y (f) 



(3.1) 
(3.2) 

(3.3) 
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Consequently, (2.9) now becomes 

(~l) m+n / (C)n(C) m f f y (-^(-m^C + z + j + l-fc) , g.fc' 
/?r(C) V n!m! (C)i(C)^!j! 1 9> 

(3.4) 

where we have used the identity (— ( — i — j)k = fc+t+j+i) k ■ The double sum of the first term on 
the right-hand side of (3.4) can be written in terms of the confluent hypergeometric function (2.6). 
In this case, we have 



< ♦ftl^WS > - (-i)-"^#S%r(- C , V. ( T 

1 + gr z <p V n\m\ g \ C, 



m 



(-i) m+ " / (0«(0m v^fv (-^(- m ) J r (C + ^ + .? + i-fc) f /y 

(3T(0 V n!m! ^ ^ ^ (CMC)^j! 1 9> 

(3.5) 

The advantage of using such expressions as (3.4) or (3.5) is that they allow us to compute the incom- 
plete gamma function only once, since its parameters are fixed. In our numerical computations of the 
variational upper bounds, we shall use (3.4) which can be implemented using symbolic mathematical 
software such as Mathematica or Maple. In order to verify the usefulness of this variational approach, 
as well as its generality for computing upper bounds for arbitrary values of the parameters A and 
g > , we will now compare our results with the exact eigenvalues obtain for certain relations of 
the parameters (as summarized in the appendix) and with the most recent approximate eigenvalues 
obtained in the literature by means of a variety of techniques. 

3.1 Comparison with exact solutions 

We first compare the variational upper bounds with the exact solutions developed in the appendix. It 
is known [24] that such exact solutions are only possible for A < and for certain algebraic relations 
between the parameters A and g . For the Schrodinger equation (2.1) with 7 = I + ^(N — 3) or 
real number 7 > — 1 and B = 1 , we have using (A. 4) for n = 1 that 

u, d ^ 1 ( 2 , 7(7 + 1) 4. 9 + 2. 9 2 (3 + 2 7 ) ^ 

H ^> = —&2 +[ r + — YTgV 2 f = ( 3 + 2 7)(1 - 25)^ (3-6) 

with the exact nodeless ground state ip{r) = r 7+1 (l + gr 2 )e~ r ' 2 l 2 . In this case, the variational 
method for computing the upper bound of H in (3.6) converges rapidly for arbitrary values of 
g > and minimization of 2 x 2 matrix (2.15) over A is sufficient to obtain very accurate upper 
bounds, as indicated by the Table 1, in comparison with the exact eigenvalues given by (3.6). 
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Table 1. Comparison of the exact eigenvalues obtained by means of (3.6) and that obtain by the 
variational technique E var discussed in the present work. The number D in the square brackets 
refers the size of the basis set used in this case. 



a 


7 


A 


-^exact 




0.1 





-0.46 


2.400 000 000 000 


2.400 000 000 000 [2] 




2 


-0.54 


5.600 000 000 000 


5.600 000 000 000 [2] 


i 





-10.0 


-3.000 000 000 000 


-3.000 000 000 000 [2] 




2 


-18.0 


-7.000 000 000 000 


-7.000 000 000 000 [2] 


10 





-640 


-57.000 000 000 000 


-57.000 000 000 000 [2] 




2 


-1440 


-133.000 000 000 000 


-133.000 000 000 000 [2] 



Further, from (A. 4) for n = 2 , we note that for A = (E± — 27 — 11) g where 

E± = 5 + 2 7 - ,g(13 + 67) ± vV(7 + 2 7 ) 2 + g(-4 + 87) + 4, (3.7) 
the correspondent eigenstates reads 

V>±(r) ^r^ +1 (l+gr 2 )(l + a ± r 2 )e- r2 / 2 , (3.8) 

where 

£ ± + (2g-l)(3 + 2 7 ) 
a± ~ (6 + 47) 

For E + , we have a + < and the wave function tp^ ( r ) nas one node on [0, 00) , thus i s 
the eigenstate for the first excited state in this case. On the other hand, for E- , a_ > and the 
wave function ^ (r) has no node and thus represent the ground state wave function when g and 
A satisfy the correspondence condition. In Table 2, we have compared our variational results with 
those obtained by (3.7). It should be mentioned that these variational eigenvalues are obtained by 
optimizing (2.15) over a iwo-dimensional subspace. The best possible value of the free minimization 
parameter A is found to be zero. 



Table 2. Comparison of the eigenvalues obtained by means of (3.7) and those obtained by the 
variational technique E var discussed in the present work. The number D in the square brackets 
refers the size of the basis set used in this case. 



A 



E ex 



0.1 -26 -15.000 000 000 000 

-12 - 1.000 000 000 000 

1 -.56174575578973345189 7.382542442102665 

-1.0182542442102665481 2.817457557897335 



-14.9999999999999 [2] 
-0.99999999999996 [2] 
7.382542442102671 [2] 
2.817457557897341 [2] 
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3.2 Comparison with the Quadrature Discretization Method 

In Table 3, we compare our upper bounds with that of Chen and Shizgal [57,58] obtained by means 
of the Quadrature Discretization Method (QDM). The QDM generally employs a discretized ver- 
sion of the Hamiltonian with respect to a set of points that correspond to the quadrature points 
associated with the chosen weight function. The distribution of the grid points is determined by 
the weight function, which controls the convergence of the eigenvalues and eigenfunctions. In Table 
3, we have compared our results with those of [57,58] using the weight functions exp(— ar 2 ) and 
exp(— r 2 i/l + A/(l + 0.5p)) , where a is chosen by the authors for the fastest convergence. The 
underlined portion of each eigenvalue indicates the convergence to the number of significant figures 
obtained by [57,58] and has been mentioned for the purpose of comparison, we may extend the 
convergence of the upper bounds obtained by the variational method by increasing the size of D 
in (2.15). The number D indicates the size of the digonalized matrix used in the diagonalization 
process or/and the number of quadrature points used in the application of the QDM. The table 
clearly indicate the faster convergence the eigenvalues by means of the variational method. 
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Table 3. Comparison of the convergence of the lowest two odd eigenvalues of the nonpolynomial 
oscillator Hamiltonian —d 2 /dr 2 + r 2 + r 2 /(l + gr 2 ) between the variational method proposed in the 
present paper and the quadrature discretization method (QDM) [57] with different weight functions 
w(r) . The number D is the square brackets refers the size of the basis set used in each method. 



9 n 


E 






w(r) — e~ ar 


2 




-Eqdm 

w(r) = e -''V 1 + A /( 1 + - 5 9) 


1 1 


3.51158094 


[D = 


1] 


3.51099389 


[D = 


10] 


3.50666367 


\D = 101 

L 1 




3.50740682 


[D = 


5] 


3.50738872 


[D = 


20] 


3 50737573 


[D = 20] 




3.50738865 


[D = 


10] 


3.50738837 


[D = 


25] 


3.50738781 


[D = 30] 




3.50738836 


[D = 


15] 


3.50738835 


[D = 


30] 


3.50738831 


[D = 40] 




3.50738835 


\N = 


16] 


3.50738835 


\D = 


301 


3.50738835 


\D = 501 


2 


7.65220193 


[D = 


1] 


7.83801690 


[D = 


10] 


7.64562064 


[D = 101 

L 1 




7.64828294 


[D = 


5] 


7.64836479 


[D = 


20] 


7 64815336 


[D = 20] 




7.64820250 


[D = 


10] 


7.64820406 


[D = 


25] 


7.64819920 


[D = 30] 




7.64820129 


[D = 


15] 


7.64820124 


[D = 


30] 


7.64820110 


[D = 40] 




7.64820124 


\D = 


19] 


7.64820124 


\D = 


35] 


7.64820124 


\D = 601 


10 1 


3.08809337 


[D = 


1] 


3.54168906 


[D = 


10] 


3 08692234 


[D = 10] 




3.08809139 


[D = 


5] 


3.08883073 


[D = 


30] 


3.08794408 


[D = 30] 




3.08809096 


[D = 


10] 


3.08809133 


[D = 


50] 


3.08805698 


[D = 50] 




3.08809088 


[D = 


15] 


3.08809085 


[D = 


60] 


3.08809057 


[D = 150] 




3.08809086 


\D = 


19] 


3.08809085 


\D = 


701 


3.08809075 


\D = 1801 


2 


7.09037623 


[D = 


1] 


11.03978681 


[D- 


= 10] 


7.08838586 


[D = 10] 




7.09037144 


[D = 


5] 


7.13612019 


[D = 


30] 


7.09012096 


[D = 20] 




7.09037060 


[D = 


10] 


7.09048160 


[D = 


50] 


7.09031284 


[D = 50] 




7.09037046 


[D = 


15] 


7.09037053 


[D = 


70] 


7.09036993 


[D = 150] 




7.09037043 


\D = 


19] 


7.09037041 


\D = 


801 


7.09037025 


\D = 1801 


100 1 


3.009831772 


[D- 


= 1] 


6.61289499 


[D = 


10] 


3.00981139 


[D = 10] 




3.009831771 


[D- 


= 5] 


3.00987806 


[D = 


100] 


3.00982785 


[D = 100] 




3.009831771 


[D- 


= 10] 


3.00983181 


[D = 


150] 


3.00983006 


[D = 170] 




3.009831771 


\D- 


= 15] 


3.00983177 


\D = 


1701 


3.00983023 


\D = 1801 


2 


7.00984496 [D = 


1] 


25.73416848 


[D- 


= 10] 


7.00981397 


[D = 10] 




7.00984495 


[D = 


5] 


7.01481472 


[D = 


30] 


7.00983900 


[D = 100] 




7.00984495 


[D = 


10] 


7.00985617 


[D = 


50] 


7.00984172 


[D = 150] 




7.00984495 


\D = 


15] 


7.00984517 


\D = 


1801 


7.00984261 


\D = 1801 
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3.3 Comparison with the Eigenvalue Moment Method 

The eigenvalue moment method (EMM) is an effective and simple technique for generating con- 
verging lower and upper bounds to the low-lying discrete spectrum of multidimensional quantum 
Hamiltonians [19,44]. The method, however, faced some difficulties in determining eigenenergies for 
certain Hamiltonian parameter values. For the non-polynomial oscillator potential (1.1), we have 
compared our variational results obtained by means of (3.15) with those of EMM [44]. As seen from 
the table, EMM is not adequate for large values of g , whereas the variational method developed 
here gives good results. 
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Table 4. Comparison of the variational upper bounds E var with the upper and lower bounds 
Eemm obatin by means of the the eigenvalue moment method (EMM) [44] . The number D in the 
square brackets refers the size of the basis set used in each method. The size of the basis set used in 
the variational method is set to be D — 18 . 



A g I E var Eemm 



0.1 0.1 1 5.186373002931507 5.1863730029314 < E < 5.1863730029316 

2 7.243961840421887 7.2439618404138 < E < 7.2439618404260 

3 9.294359110874627 9.29435911086337 < E < 9.29435911088159 



0.1 0.5 1 5.100857624300696 5.100842 < E < 5.100865 

2 7.11898087156427 7.11890 < E < 7.118901 

3 9.131812401691521 9.131799 < E < 9.131838 



0.1 1 1 5.065569521783354 5.06428 < E < 5.06609 

2 7.073726361909647 7.0730 < E < 7.0744 

3 9.078911720303639 9.0787 < E < 9.07892 



0.5 1 1 5.893595152339402 5.89359515233919 < E < 5.89359515233945 

2 8.1778716934677 8.177871693435 < E < 8.177871693485 

3 10.429204118147453 10.4292041181366 < E < 10.4292041181548 



1 0.1 1 6.704238892478644 6.7042388924777 < E < 6.7042388924788 

2 9.261914780826569 9.2619147807 < E < 9.2619147809 

3 11.760620962669972 11.7606209626312 < E < 11.7606209626917 



1 1 1 5.65139331725017 5.6503 < E < 5.6521 

2 7.73482804292358 7.734 < E < 7.736 

3 9.787669778509466 9.7875 < E < 9.7881 



100 100 



1 

2 
3 



5.993438873366758 
7.996024673021835 
9.997153638602487 



- < E < 6.389 
7.9947 < E < 8.037800 
9.9969 < E < 10.0113 
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3.4 Some remarks about first-order perturbation expansions 

The closed-form expressions for the matrix-elements (2.9) can be used explicitly to express the first- 
order perturbation expansions for the perturbed Hamiltonian (2.1) for any dimension N . Indeed, 
we have by means of (2.9) that 

In particular, if A = and S = 1 , we have 

In Table 5, we have expressed the first few states for different values of g and for small values of A 
in the case of N = 3 and 1 = 0. 

Table 5. The first-order perturbation expansions of the lowest three eigenvalues for different values 
g and valid for small values of A . For comparison see Ref. [28] . 

JL ^oo - 3 Ef Q - 7 ff|o -11 

0.5 0.741 907 668 608 869A 1.058 912 626 973 532A 1.216 304 645 477 702A 

1 0.515 744 312 282 624A 0.648 934 634 510 945A 0.712 059 158 182 293A 

2 0.327 839 771 209 399A 0.374 239 389 891 732A 0.396 830 711 146 285A 
5 0.160 824 698 125 120A 0.169 855 585 054 398A 0.174 677 060 833 675A 
10 0.088 111 301 584 084A 0.090 376 621 370 169A 0.091 681 154 216 976A 
20 0.046 566 260 901 426A 0.047 083 627 877 117A 0.047 401 822 090 664A 
100 0.009 831 778 572 526A 0.009 844 958 882 145A 0.009 853 945 311 642A 
500 0.001 992 603 359 200A 0.001 992 880 765 846A 0.001 993 079 888 610A 



4. Conclusion 

We have developed and applied a variational method to compute accurate upper bounds for the 
energies corresponding to the non-polynomial oscillator potential (1.1) in arbitrary dimension N. 
The results show the effectiveness of this approach in comparison with other techniques. We have 
compared our variational results with those of many other techniques developed in the literature, 
for example with the Pseudo-Perturbation Expansion Method [53], the Shifted l/N Expansion [47], 
and the matrix continued fraction algorithm [49] . For each technique, we have confirmed their reults 
with the use of very few matrix elements. Indeed, all of our variational bounds were obtained by the 
optimization of a square matrix (2.15) with size D not exceeding 20 ; these results can, of course, 
be further improved by simply increasing D . 
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Appendix: An infinite set of exact solutions 



In this appendix we investigate the exact solutions of the non-polynomial oscillator Schrodinger 
equation in the N -dimenional case. The existence of an infinite number of exact solutions of the 
one-dimensional Schrodinger equation (2.1) (B = 1,N = 1,1 = 0) with eigenfunctions given by 
products of exponential and polynomial functions of r 2 , for specific relations between the couplings 
g and A , was first reported by Flessas [22] and soon afterwards was extended by Varma [23] . 
The main point of exploring these exact solutions in this appendix is that: (1) they generalize 
the previous exact solutions obtained for the one-dimensional case and they are valid for arbitrary 
7 = I + i(7V — 3) > — 1 ; (2) they provide interesting examples for testing the vadility and accuracy 
of our variational technique. Assume 



,.2'i 



and substitute it into equation (2.1) yields ( B = 1 , E = E^ ) 

{E - 2 7 - 3)a + (6 + 4 7 )ai + ^ j ( - 2g(2i - 2) - A + g(E - 2 7 - 3)j a,_i 

+ [2ig{2i - 2) - 4i + 4i.g( 7 + 1) + E - 2 7 - 3 ] a, + [2(i + l)(2i + 2 7 + 3) ] a l+1 Xr 21 



(Al) 



+ | ^ - Ag{n - 1) - A + g(E - 2 7 - 3) J a„_i + \ [ 2ng(2n - 1) - 4n + 4n 9 ( 7 + 1) + E - 2 7 - 3j a n 
+ (g(E - 2 7 - 3) - A - 4 5 n) a n r 2n+2 = 

The last term of this equation implies 

A = #(£-4n-2 7 -3), n = 0,1,2,... (A2) 

provided that the energy E and the interaction parameters A, g satisfies the (n + 1) x (n + 1) - 
tridiagonal determinant 
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where 

( a n = g(E - 4n - 2 7 - 3) - A, 

b n = 2ng(2n + 2 7 + 1) - An + E - 2 7 - 3 = E - 4n - 2 7 - 3 

v c - 2n(2n + 2 7 + 1), rc = 0, 1, 2, . . . 
This detcrminantal equation can be written using the LU-decomposition as 
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in which 



c n a n —\ 

u n = o n , u = b , v n = a n , t n 



U n -1 Un-l 

This immediately leads to the following simple expression of computing the polynomial conditions 
as 

n n / > 

^-n^n— 1 



ti n / 

Vn(E) - l[u k = H (b„ 
fe=l k=l ^ 



U n -1 



0. 



(AA) 



It should be clear that our relations (A.3)-(A.4) summarize the exact solutions not only in the 
N -space dimensional where 7 = I + ^(N — 3) but also for arbitrary real values of 7 > — 1 . 
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